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ABSTRACT 
The  general  dissipative  Braginskii  single«fluid  model  is  applied 
to  simulate  tokamak  transport.  An  expansion  with  respect  to 
e  ■  («i>«Tj)  ,  the  factor  by  which  perpendicular  and  parallel  transport 
coefficients  differ,  yields  a  numericlly  tractable  scheme.  The 
resulting  1  1 /2  D  procedure  requires  computation  of  2D  toroidal 
equilibria  with  flow  together  with  the  solution  of  a  system  of  ordinary 
1D  flux-averaged  equations  for  the  time  evolution  of  the  profiles. 
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I.   INTRODUCTION 

As  tokamak  experiments  become  larger,  cleaner  and  hotter, 
realistic  simulation  of  particle  and  energy  transport  is  becoming 
essential.  Such  a  model  is  developed  here  in  the  context  of 
single-fluid  theory  in  conjunction  with  numerical  techniques  suitable 
for  representing  neoclassical  tokamak  transport.  The  equations  of 
dissipative  magneto-hydrodynamics  with  axisymmetry  imposed  and  with 
source  terms  are  employed.  Tensor  viscous  effects  are  included  in  the 
equations  of  momentum  and  energy  conservation,  anisotropic  heat  flow 
and  electrical  conductivity  in  the  energy  equation,  and  anisotropic 
electrical  conductivity  and  Hall  effect  in  Ohm's  law.  The  general 
Braginskii  forms,  Ref  [1],  of  the  viscous  stress  tensor,  anisotropic 
heat  flow  vector,  and  anisotropic  conductivity  tensor  are  used. 
Neutral^bearmheated  experiments  have  shown  large  toroidal  plasma  flow 
comparable  with  the  ion  sound  speed  but  no  significant  induced  poloidal 
flow;  this  damping  of  the  poloidal  flow  -*  even  in  the  case  of 
perpendicular  injection  -*  is  due  to  enhanced  transport.  Our  model 
aims  to  study  the  effect  of  large  toroidal  flow  on  transport  due  to 
viscosity,  which  tends  to  couple  this  large  flow  to  poloidal  flow. 
This  inclusion  of  flow,  together  with  the  Hall  term  makes  our  model 
substantially  different  from  other,  simpler  neoclassical  models. 

A  rigorous  treatment  is  applied  to  this  set  of  general 
single-fluid  equations  using  classicalMtype  transport  coefficients  and 
its  properties  are  discussed.  This  framework  may  be  inadquate  to 
simulate  anomalous  transport,  but  we  definitely  prefer  an  exact 
treatment  to  a  purely  phenomenological  model.  This  is  true  especially 
because  a  generalization  of  the  pressure  tensor  beyond  the  Braginskii 
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form  is  not  yet  known.  At  the  same  time  when  empirical  approaches 
aimed  at  the  interpretation  and  prediction  of  experimental  results 
become  more  complicated  by  the  inclusion  of  toroidal  geometry,  it  is 
very  desirable  to  understand  the  properties  of  the  classical  model. 
Therefore,  many  of  the  relations  derived  here  should  be  of  general 
interest  and  should  help  to  improve  empirical  models. 

In  a  previous  paper  [2]  we  discussed  the  properties  of  this 
general  fluid  model,  especially  with  respect  to  numerical  solution.  A 
major  difficulty  arises  because  the  transport  coefficients  are  grossly 
disparate  in  magnitude,  so  that  for  typical  tokamak  parameters  the  very 
small  perpendicular  flow  most  likely  cannot  be  computed  from  the 
general  equations.  We  now  solve  the  equations  for  the  different  scales 
step  by  step  by  applying  an  expansion  in  e  =>  (<djt^)  ,  where  w*  is  the 
gyrofrequency  and  t^  the  collison  time  for  ions.  The  leading»-order 
transport  coefficients  are  the  heat  conductivity  and  viscosity  along 
the  magnetic  field,  which  cause  the  temperature  to  be  constant  along 
the  field  and  allow  only  toroidal  flow.  The  equilibrium  is  given  by  a 
generalized  Grad^Schltiter'-Shafranov  equation.  The  next-=order 
correction  yields  a  flow  along  the  magnetic  field  and  then  the 
second-order  correction  determines  diffusion.  From  this  second-order 
solution  the  flow  perpendicular  to  the  magnetic  surface  is  the  main 
interest.  Hence  a  significant  simplification  is  achieved  if  the 
surface  average  of  the  equations  is  taken,  a  procedure  which  determines 
the  diffusion  entirely.  In  this  fashion  self-consistent  steady^state 
solutions  are  obtained  where  the  outflow  is  governed  by  the  equilibrium 
field,  the  sources  and  the  transport  coefficients.  The  fact  that  large 
toroidal  flow  is  consistent  with  viscosity  and  its  influence  on 
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diffusion  have  alredy  been  discussed  in  Ref.  [33.  However,  the  model 
used  there  neglects  Ohm's  law  and  hence  the  corrections  on  the  fields 
and  is  therefore  less  complete  than  ours. 

To  get  a  more  complete  picture,  profile  changes  on  the  diffusion 
time  scale  are  included,  so  that  all  profiles,  are  fully  determined  by 
sources  and  diffusion.  We  follow  the  procedure  given  by  Pao  in 
Ref.  [*G  to  determine  the  diffusion.  The  resulting 
integro-differential  equations  for  the  profile  functions  are  quite 
complicated,  but  can  be  simplified  by  expressing  the  time  evolution  of 
the  profiles  in  terms  of  the  evolution  of  the  flux  function.  Then  only 
that  flux  function  determined  by  the  time  derivative  of  the 
Grad-SchluterrShafranov  equation  needs  to  be  determined  as  a  function 
of  r  and  z;  the  remaining  dependence  is  one'-dimensional.  This  scheme 
alternating  between  2D  equilibria  and  1D  transport  is  quite  similar  to 
Grad's  1  1 /2  D  algorithm  [5-7].  However  our  model  is  made 
significantly  more  complicated  by  the  inclusion  of  flow  and  the  more 
general  dissipative  fluid  model.  Our  ansatz  for  simulating 
two-dimensional  tokamak  transport  is  similar  to  that  of  Hirshman  and 
Jardin,  Ref.  [8],  but  is  made  more  general  by  the  inclusion  of  the 
full  viscous  tensor  together  with  flow  in  the  basic  equilibrium,  which 
leads  to  a  poloidal  variation  of  the  pressure  and  density  on  magnetic 
surfaces.  Also  the  application  of  the  f inite*element  method  to  the 
solution  of  the  equilibrium  equations  yields  a  completely  different 
numerical  scheme. 

In  our  approach  to  the  simulation  of  the  long-time  behavior  of 
tokamaks  the  transport  coefficients  are  required  for  the  macroscopic 
description.  The  transport  coefficients  are  assumed  to  be  given.   In 


-5- 
the  collision-dominated  regime  they  are  explicitly  known.  The  various 
transport  coefficients  are  of  greatly  disparate  magnitude.  The 
parameter  e  =  (oojij)^1  measures  the  relative  smallness  of  different 
elements,  e.g.,  the  heat  conductivity  along  the  magnetic  field  B, 
compared  with  that  in  poloidal  direction  but  perpendicular  to  B,  and  is 
in  the  range  of  10*"3  to  10  ,  being  10~5  for  hot  tokamaks.  Such  values 
are  so  small  compared  with  the  inverse  aspect  ratio,  plasma  beta  and 
B  /Bt  that  we  take  all  the  latter  quantities  of  order  one  in  e,  and  we 
may  then  take  the  neoclassical  corrections  to  Braginskii  transport  all 
to  be  factors  of  order  one.  In  the  neoclassical  regime  the  trapping  of 
the  particles  is  taken  into  account.  The  multiplicative  correction 
factors  to  the  standard  terms  for  the  representation  of  neoclassical 
transport  yields  local  transport  coefficients  which  may  require 
non-local  averaging  for  their  determination  and  may  depend  on  the 
equilibrium.  The  simulation  of  plasmas  by  one-dimensional,  transport 
codes  and  direct  experimental  observation  shows  that  some  of  the 
transport  coefficient  are  anomalous  -ft  possibly  owing  to  instabilities, 
fluctuations  or  turbulence.  The  perpendicular  electron  heat 
conductivity  and  the  diffusion  factor  exhibit  anomalous  behavior;  their 
values  exceed  the  classical  terms  by  a  factor  10  to  10J  as  shown  in 
Refs.  [9,10].  It  is  emphasized  that  such  a  correction  factor  of  the 
order  of  e  is  almost  tolerable  in  our  systematic  expansion 
procedure.  A  more  phenomenlogical  choice  of  the  transport  coefficients 
in  certain  limits  is  thus  possible,  especially  if  the  value'  of  o>t  is 
chosen  freely  to  match  experimental  data  with  the  sole  restriction  that 
it  be  large  enough  to  remain  consistent  with  the  expansion.  In  our 
opinion  it  is  more  important  to  examine  the  two-dimensional  properties 
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of  the  equilibria,  especially  the  toroidal  flow,  as  well  as  the  effects 
of  viscosity  and  the  Hall  constant  in  the  first  place  and  then  to  adopt 
the  magnitude  of  anomalous  factors  for  agreement  with  experiments. 

The  plan  of  the  paper  is  as  follows:  Section  II  presents  the  fluid 
equations  and  introduces  the  basic  quantities  needed  for  the  analysis. 
Then  a  systematic  expansion  in  e,  which  measures  the  magnitude  of 
different  transport  terms  and  defines  the  natural  scale,  is  employed. 
Section  III  gives  in  lowest  order  in  e  the  steady^state  solution  in 
terms  of  a  Grad-SchlUter-Shafranov  equation  for  the  poloidal  flux 
function  but  with  the  pressure  as  a  specified  function  of  r  and  z,  and 
a  purely  toroidal  flow  whose  angular  velocity  is  a  function  of  the  flux 
alone,  as  is  the  temperature.  In  Section  IV  it  is  found  that  in  first 
order  in  e  a  poloidal  flow  is  required  to  achieve  momentum  balance,  but 
no  particle  or  energy  outflow  occurs.  In  second  order  in  e,  presented 
in  Section  V,  particle,  momentum,  energy  and  current  sources  are 
admitted.  The  time  evolution  of  lowest-order  equilibrium  profiles  is 
also  included.  The  equations  are  only  solved  to  the  extent  necessary 
for  determining  the  slow  temporal  evolution  of  the  lowest-order 
equilibrium  quantities.  Consequently,  flux-surface-averaged  equations 
are  obtained  for  these  quantities.  In  particular,  the  energy  flux  and 
the  evolution  of  the  temperature  are  substantially  different  from  the 
simpler  disslpative  MHD  models.  The  discussion  and  conclusions  are 
presented  in  Section  VI. 

This  set  of  self-consistent  equations  for  describing  the  temporal 
evolution  of  tokamaks  on  the  slow  diffusive  time  scale  is  presented  in 
a  form  suitable  for  numerical  treatment.  The  basic  2D  equilibrium  and 
the  Grad-Schluter^Shafranov  type  equations  required  later  may  be  solved 
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by  finite-element  methods.  The  code  of  Kerner  and  Jandl,  Ref.  [11], 
yields  the  equilibrium  in  flux  coordinates,  thus  making  computation  of 
the  required  surface  integrals  straightforward.  The  discretization  of 
the  time  evolution  leads  to  systems  of  ordinary  diferential  equations. 
Our  analysis  includes  the  basic  concepts  for  accurate  and  efficient 
numerical  solution.  Details  of  the  numerical  method  and  the  code 
itself  will  be  described  elsewhere. 
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II.   PHYSICAL  MODEL 

The  plasma  is  described  by  a  macroscopic  single-fluid  model.   This 

model  relates  the  mass  density  p,  velocity  v,  scalar  pressure  p, 
pressure  tensor  t,  specific  internal  energy  e  (i.e.,  energy  per  unit 
mass),  and  the  magnetic  field  B: 


Continuity:  —  p  +  V-(pv)  =  ps  (1) 

at   -  - 


Momentum:    p(—  +  v«v)v  =  >-   Vp  +  JxB  +  V«ff  +  Ps  (2) 

9t  = 


Energy:      p(—  +  v«V)e  +  pV«v  =  -  V«£  +  I    :V  v  +  J«(E  +  vxB)  +  Es  ,   (3) 
9t  = 


where  ^  denotes  the  heat  flux,  and  ps,  P_s  and  E3  are  sources  of  mass, 
momentum  and  energy.  The  heat  flux  is  of  the  form 


K/\ 
£  =    K-VT  =    KV      T    +    KVT   -   BxVT    ,  (4) 

II     II  11  B 


K  being  the  thermal  conductivity  tensor  and,  in  particular,  K   the 

I  I 
parallel  and  K/\  and  K  the  perpendicular  conductivities.   The   ideal 

gas  law  assumes  the  form 


Ideal  Gas:  p  =  pRT  ,  (5) 

where  the  constant  R  is  the  gas  constant  RQ  divided  by  the  effective 
molecular  weight  of  the  gas. 
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We  utilize  the  first  law  of  thermodynamics 

TdS  =>  de  +  p  dx  =  de  +  p  d(1/p)  (6) 

to  replace  the  specific  internal  energy  by  the  specific  entropy  S  or 
the  temperature  T  and  obtain  an  equation  for  the  specific  entropy: 


p  T(_  +  v«V)S  =»  -V«2  +  ff:V  v  +  J«(E  +  vxB)  +  Es  -  psRT  .     (3') 
3t  = 


For  an  ideal  polytropic  gas  the  specific  internal  energy  can  be 
expressed  in  terms  of  the  temperature  and  the  specific  heat  at  constant 
volume  as: 


e  =  cvT  ,  (7) 


and  the  energy  equation   in  terms  of  T  reads 


<v3 


p(— :  +   v*V)T  +   pV-v  =  -   V.q  +   I    :V   v  +   J«(E   +   VxB)   +   E3    .      (3") 

ot  =  —     —       —  — 


The  entropy  for  an  ideal  polytropic  gas  is  given  up  to  an  arbitrary 
constant  as  [12] 

S  -  SQ  =  R  an(T1/Y"1/p)  =  cy  Sln(Pp"Y)  ,  (8) 

where  Y  is  the  ratio  of  specific  heats.   The  relation  between  the 
electric  field  E  and  the  current  J  is  given  by 
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s 


Ohm's  law:       n«J  +  nHJ.«B  =  E  +  vxB  -  Jc  , 

with  the  resistivity  tensor  n,   the  Hall  constant  nH  and  a  current 
source  J3.   It  is  convenient  to  write  the  mass  density  explicitly  for 

the  Hall  term  as  n^  =  ^h^*3  =  —  /'p* 

Maxwell's  equations  close  the  system: 

3B 

—  =  -  VxE  , 
Maxwell:  3t     

J  =  VxB  ,  (10) 

V'B  =  0  . 

In  Ref.   [2],  eq.   (63)  we  established  the  following  form  for  the 
stress  tensor: 


B  B 
yQW  +  u, (W-B  -  B.W)  +  U2(B.W.B  -  B'W-B) 


B2 


B  B 
+  uJlI  -  3  JB-W-B  +  un    (B  B.W'B  -  B^W-B  B)  ,   (11  ) 


and  the  condition  for  the  viscous  entropy  to  be  positive  (Ref.   [2], 
eq.   (64)): 


U2  +  3u3  +  B2y5  =  0  .  (12) 


It  was  also  shown  that  this  form  agrees  with  that  of  Braginskii  if  the 
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viscous  coefficients  u  (v  =  0,1,2,3  and  H)   are  related  to  those  of 
Braginskii  by 


u0  =  n2  . 

»2   -  Cn-i  *   n2)/B2  ,  (13) 


u3  --  (n0  +  n-|  -  2n2)/2B2 
u4  -  (2nu  -  n3)/2B3  . 


The  viscous  stress  is  related  to  the  strain  matrix  S: 


3v,   3v, 


V-v  =  — i  +  — ±  -  4  6n  v'v  •      (M) 

J 3x,    3x,    3   1J 


=ij   ,ij   3  iJ  -  -   3Xj   3Xl 


The  field  tensor  B  satisfies  the  relation  with  any  vector  a: 

B^a  -  Bxa  and  a«B  =  axB  .  (15) 

The  usual  cylindricl  coordinates  r,9,z  are  used.   The  velocity  is 
represented  by  its  contravariant  components: 

v  -  uV  r  +  vV  6  +  wV  z  .  (16) 

For  axisymmetric  flows,  3/36  ■  0,  we  have 


V-v  -  1  [(ru),n  +  (rw)  „]  .  (17) 
r 
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where  ,r  denotes  the  partial  derivative  with  respect  to  r,  etc. 
In  vector  notation  the  strain  matrix  assume  the  form 


S  =  V  v  +  (V  v)tr  =  VuVr  +  VrVu  +  r[v(I)Ve  +  VeV(-j] 

=  r  ~         ~  ~  r 


-  I   [VrVe  +  VeVr]  +  2urV0Ve  +  VwVz  +  VzVw  .  (18) 

r 


The  magentic  field  is  divergence-free.  It  is  convenient  for 
axisymmetric  fields  to  introduce  the  poloidal  flux  \p  and  the  poloidal 
current  profile  x: 


b  =  v^xVe  +  x  ve  •  (19) 


The  current  then  becomes 


J  =  7X  x  V6  -  J  S79  ,  (20) 


where 


Jn  =  r2V«(Vi|>/r2)  =■  A%  .  (21  ) 


o 
The  Lorentz  force  is  then  given  by 

JxB  =  'JQ/r2  7ij;  *•  x/r2  Vx  +  V6VxxViJ;«?6  .  (22) 

As  stated  above  the  various  transport  coefficients  are  of  greatly 

disparate  magnitude  (see  Ref.   [11]).  The  parameter  e  =  (oo_  •t<)'" 

ci  1 

measures  the  relative  smallness  and  is  in  the  range  of  10"'  to  10   , 
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being  e  r  10  for  typical  tokamak  parameters.  Such  values  are  so 
small  compared  with  other  plasma  parameters  such  as  the  inverse  aspect 
ratio  and  beta  as  well  as  the  neoclassical  corrections  to  Braginskii 
transport  that  we  take  all  these  to  be  factors  of  order  one.  The 
scaling  of  the  transport  quantities  is  therefore  established  by  the 
classical  expressions.   From  the  definitions  it  follows  direct  that 


K    :  KA  :  K  -  1  :  £  i  £2 


and 


2    2 
no  :  n3  :  n4  :  n1  :  n2  =  1  :  e  :  e  :  e 


£ 


From  Ohm's  law  it  is  found  that  the  parallel  or  perpendicular 
resistivity  are  a  factor  of  about  e  smaller  than  the  Hall  term. 
Comparing  the  dissipation  terms  in  the  energy  equation  reveals  that 
perpendicular  resistivity  and  parallel  heat  conductivity  are  both  of 
order  e  .  The  following  ordering  is  then  given: 


K   ,u3  -9(1)  . 

nH'KA1  »u1  ,y4  =  e^  •  ^23) 


K  ,n  ,n. .u0,w2  =  9(e  )  . 

We  now  employ  a  systematic  expansion  in  e.   To  describe  diffusion,   all 
quantities  are  expanded  up  to  second  order  f  =  f^°'  +  zf         +   e  f   . 

The  general  fluid  equations  in  Braginskii  contain  a  friction  force 
in  the  momentum  equations  for  ions  and  electrons,   which  cancel, 
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however,  if  both  equations  are  added,  and  a  heat  generation  term  in  the 
energy  equation  as  a  consequence  of  collisions  of  electrons  with  ions 


,  ,   ,      BxVT 


where  equal  temperatures  Tg  =  T^  is  assumed.   The  comparison  of  Qe  with 

p 
the  Ohmic  dissipation  n  J~"  reveals  that  both  terms  up  to  a  factor  of 

beta  are  of  the  same  order.   Consequently  Qe  contributes  only  in  second 

order. 
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III.   ZERO  ORDER 

In  leading  order  in  e  there  are  no  source  terms  and  we  have  to 
look  for  steady-state  solutions  3/3t  =  0. 

First  the  energy  equation  is  discussed.   The  basic  equilibrium 
must  conserve  the  entropy.   From  eq.   (3')  we  have 


;(o) 
dt~ 


p(o)  T(o)  dS -o-V-CkI^V   T(o)]  ♦  (I,  V  v)(0)  .    (24) 


The  pressure  tensor  reads 


B(o)B(o) 

!(0)  =  u{o)[U   -  = = )  B(o).W(o).B(o) 

3   V      B(o)2  '  ~ 


The  energy  equation  then  yields 

V   T  =  B(o)-VT(o)  -  0 

which,  because  of  axisymmetry,  implies  that 


.(o)  .  T(o)(|jJ(o))  %  (25) 


The  viscous  entropy  production  is  given  by 


(E:V  v)(o)  =  1  ff(o).W(o)  =  M^o)(B(o)-W(o).B(o)]2  , 


with  uio)   =  -3uio)/B(o)2,  and  hence 


5  '2 
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,(o)  .  w(o)  .  B(o)  .  0  (  (26) 


with  the  non-trivial  solution 


v(o)  =  r2jj(o)u(o))Ve  f  (27) 


since  then 

w(o)  .  r2n(o)'(y^(o)ye  +  veV*(o))  ,        (28) 

corresponding  to  zero  parallel  viscous  stress. 

The  prime  denotes  differentiation  with  respect  to  the  argument. 
The  continuity  equation  is  trivially  satisfied.  Ohm's  law  and 
Faraday's  law  are  satisfied  by  introducing  a  potential 


;(o)  .  _70(o)  .  _v(o)xB(o)  .  _a(o)v^(o)  f  (29) 


The  component  of  Ohm's  law  along  B^0'  gives 


,(o)  .  4,(0)^(0))  (30) 


and   the  component  along  7i|/0' 


*♦£!-  G(o)    .  (3D 

d*(0) 


In  this  order  the  energy  production  and  the  viscous  stress  tensor 
vanish  so  that  the  momentum  balance  reads 
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p(o)v(o).yy(o)  =  -p^o)rn(o)2Vr  =  -Vp(o)  ♦  j(o)xB(°: 


The  6-component  implies,  using  eq.   (22),  that 


„,(o)  .  „,(o)(i|>(o))  _  (32) 


The  B^0'  component  implies,  when  one  uses  the  equation  of  state, 


B<0>.V[r2a(,0)f  -  in  p(°»]  -  0  . 
2nT(o) 


The  density  is  then  given  by 


p(o)  -  ;<o)(#(o)j  exp  [d^l)    ,  (33) 

2RT(o) 


with  p^0'  an  arbitrary  constant  on  flux  surfaces. 

Note  that  both  the  density  and  the  pressure  vary  on  a  flux  surface 
for  non^zero  toroidal  flow,  i.e.,  p(o)  =  p (o) U(o) ,r) .  With  eq.  (33) 
the  momentum  equation  is  of  the  form 


3P(01  vj<o)  .  j(o)xB(o)  m  (31|) 


9*(o)  " 


The  final  equilibrium  equation  is  given  by  the  Vi/0'  component  of  the 
momemtum  blance 
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aVo)  .  -r2  1P|!1_  x(o)  (o)'  .  (35) 


This  is  the  Grad-Schluter-Shafranov  equation  for  the  flux,  where  the 
pressure  contains  a  contribution  from  the  toroidal  flow. 

We  now  collect  the  results  from  zero  order.  There  are  four  free 
surface  functions,  namely  T(o),  ft(o),  x(o)  and  p(o).  With  a  specified 
plasma  shape  the  poloidal  flux  is  evaluated  from  the  elliptic  PDE,  eq. 
(38),  with  the  boundary  condition  that  \\>^°'  be  constant  on  the  plasma 
surface.  It  is  further  required  that  the  pressure  vanish  and  the 
density  be  constant  on  this  surface.  The  geometry  of  this  magnetic 
field  configuration  is  then  used  in  the  averaging  process  to  determine 
diffusion  later  on.  All  other  quantities  are  then  defined.  The 
entropy  of  this  system  is  given  up  to  an  arbitrary  constant  according 
to  eq.   (8): 

s(o)  .  R  tnrfoH/y-ly^o))   m  (36) 

The  PDE,  eq.  (35),  is  solved  by  the  finite-element  method  where 
higher-order  elements  (8-node  quadrilateral  elements)  are  used.  By 
mesh  rearrangement  in  the  iteration  the  nodel  points  are  put  on 
surfaces  of  constant  flux.  The  final  equilibrium  is  obtained  in  flux 
coordinates  without  directly  employing  a  flux  coordinate  system.  In 
this  fashion  equilibria  are  computed  very  efficiently,  as  demonstrated 
in  Ref.  [11]. 
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IV.   FIRST  ORDER 

The  firsts-order  corrections  in  e  to  the  basic  equilibrium  are 
evaluated.  Again  there  must  be  no  source  terms  and  no  time  variation 
((3/3t)  -  o).  We  start  with  the  continuity  equation,  which  is  of  the 
form 

V-(pv)(1 }  =  0  =  V.(p(o)v(1 })  ,  (37) 

which  is  satisfied  by  the  introduction  of  a  stream  function  for  the 
first-order  velocity: 

p(0)v(,)  =  vx(1)xve  +  P(0)v(1)ve  .  (38) 

Next  Ohm's  law  is  examined,  which  assumes  the  form 


^H 


(D/p(o)  3P^  7(|)(o)  m   ,v  (1)  +  v(o)xB(1)  +  v(1)xB(o)  t        (39) 


3* 


(o)  "T 


The  6-component  yields 


ve.v<1)xB<°>  =  o  -       1       b(q).vx(1)   , 

-  -    -         r2p(o)  -    - 


which  implies  for  finite  p^ 


A(1)  =  x(1)(^,(o))  §  (J40) 


The  projection  along  B^0'  gives 
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B(o).V[*(1)    -    a<°V1>]    =    0    , 


with  the  solution 


4,(1)   =   AH*™)    *   a(oV1}    ,  (41) 


where  4>q  is  a  free  surface  function.  The  r— component  of  Ohm's  law 
determines  the  toroidal  component  of  the  velocity.  We  then  obtain  for 
the  velocity 


V(D  .  iilll  B(o)  +  ve  r2[n(D/p(o)  ar>l  +  a(o)'„0)  +  (1)']   (42) 

p(0)  ~  3iJ-(0) 


which  states  that  apart  from  a  toroidal  component  the  velocity  is 
located  within  a  flux  surface.   The  energy  equation  is 

p(o)T(o)y(1).£S(o)  =  y.(KVT)(1)  +  (ff:V  v)(1)  +  J. (E  ♦  v  *  B)(1)  . 

The  left-hand  side  is 

LHS  =  T(o)A(1),B(o).VS(o)  . 

The  first  term  on  the  right-hand  side  only  contributes  through 


K  (o)  ^  ^ 

RHS1  =  B(o)-V[   II    B(0).V(T(1)  -  T(o)V1))  +  *A    >'t(o)'] 
-  B(o)2  -    ~  B(o) 
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The  firstiorder  viscous  entropy  production  vanishes  according  to  eq. 
(26)  and  the  third  term  only  contributes  in  second  order  in  e. 
The  energy  equation  then  yields 


B(o).v{_T(o)x(1)'s(o)  +  A    (o)T(o)« 


(1) 
A 


B(o) 


K(o) 


♦   II   B(0)-V[T(1)-T(0)'*(1)]}  -  0  ,      (43) 


B(0)2  " 


which  can  be  integrated 


_1L  B<°>.V[T(1)-T<°>yi>] 
B(o)2   " 


K.(1) 


=  -      A        x(o)T(o)'    +   T(o)x(1)'s(o)    +   F(l^(o))    i      (144) 


B(o) 


where  F(iK°M  is  a  £ree   surface  function.   The  integration  constant  is 
now   determined   from   the    requirement    that    the    function 
f(U  =  T(1)  _>  t^'i^/1^  be  single-valued  and  hence  that  the  integration 
over  a  closed  flux  surface  vanish. 
The  operator 


B(o).Vf  =  V4)(o)xVe«Vf 


is  evaluated  in  a  flux  coordinate  system  i/o),9,6p  with  the  Jacobian 
V4)(o)x70-V6p  =  1/J,  which  yields 
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,(0).Vf  =  lil=  B<°>  *£ 


(1) 


J  86 


when  the  arc  length  61   on  a  surface  4>^0'   =  ct.   in  a  plane  9  =  ct.   is 
introduced.   For  nont-zero  B^0)  eq.   (44)  then  reads 


-(o) 


L  B<o)  ±   [T(1)-T(o)V1)] 


,(o)2   P   dJL 


K  <1) 
.  A    (o)T(o)*  +  T(o)x(1)'s(o)  +  F(U/(o)) 

B(o) 


The  integral  over  a  closed  magnetic  surface  is 


0  =  _  (o)T(o)' 


d£      _  B(o) 


B(o)  K(o) 


(o)2 


+  T(0)A(1)'  I   «   B^SCO)  +  r       a       b 

T  B(o)    K(o)  J  B(o)  K(o) 


This  expression  is  divided  by 
function  g  is  defined  as 


dg, 
,(o) 


and  the  surface  average  of  a 


<g> 


da. 

>(o) 


g  /  o 


dg. 

B(o) 
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Note  that  this  definition  gives  the  usual  Pf irsch-Schluter  corrects 


to  the  transport  coefficients.   The  function  F ( vi^ 


(o) 


is  determined  as 


,(o)2 


f(*<°>)  <^iz>  =  x(o)T(o)'  <; 


•  (o) 


c  (1) 
V\ 

c  (o) 


B(o)>  .  T(o)xdV  <B^s(o)> 

lf(0) 


This  result   is  inserted  in  eq.   (44),  which  is  again  integrated.   The 


first'-order  temperature  is  then  given  as  an  integral  along  \p 
from  an  arbitrary  starting  point  up  to  the  actual  point  P: 


(o) 


ct. 


T(D  =  T(o)yn  +  T(nu(o)}  +  (o)T(o)' 


x  { fp  jl.  *i!li  ^ 

JO   R(o)  ^(o)    .R 


(1)r(0)/k.(0) 


B^u'/Klu'> 


<B(o)2/K(o). 


dJ, 


v  (1) 


o  B(o)  „(o: 
Bp   kii 


Bl 


+  T(o)A(1 )' 


dl      S(0)b(°)2    fP  di      B(o>2 


0  B(o)     K-(°) 
BP       KM 


<B(o)2s(o)/K(o: 


0  Bp0)  K(o)    <b(°)2/K(0)> 


-}  ,   (45) 


with  an  arbitrary  constnt  T(1)(^(o))  on  each  surface.  The  momentum 
equation  reads  in  first  order 

pv-V  v(1 ^    +  Vp(1 )  -  (JxB)(1 ^  +  V«ff(1 ^  . 


The   left-hand  side   is 
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LHS  =  Vp<1>    ♦   Vr[^p(1)m(0)2   -   2n(°V1Vo)/r]    +   VeA(1  }  V0^  .Vr2fi(o) 


We  calculate  the  stress  tensor 


B(o)B(o) 
E  (1)  .  U(D(W(0).B(0)  -  B(o)-W(o))  +  u^0)(H  "  3  ~   ~    )  (B-W.B)(1) 


,(o)2 


+  41)(B(o)B(o).W(o).B(o)  -  B(o).W(o).B(o)B(o)) 


From  eq.   (28)  we  already  have  W(o)  and  with  eq.    (18)   we  evaluate 
W(1 K      We  then  get 


B.W.B<1>  =  .2x(o)fl(o)»B(o).Vt(1)  -  <  B(o)2  iiLil  B(o)  rj^ 
~  =  "  -  3       p(o)   r   RT(o) 


+  £^  B<0>.V[B<°>2]  +  2X(o)B(o).V  [!£_  &H  *    G(o)V1)]  .   (46) 
p(o)  "  "      Lp(o)  3^(o) 


The  part  of  the  first-order  stress  tensor  proportional  to  u^0'  is 
defined  by  ffj1  ) : 


B(o)B(o) 

ff(1)  -  fii  -  3  Z-—L )  n(1)  (47) 

'  B(o)2 


where 


->       -    =    —  -5     „(0)   _ 


+  2u^rB(o)[x(o)n(DRT(o)  [^1]  2  B(o)2  X™*   O^ij  .   (M) 

RT(o)  >(|/°;   3       p(o)  RT(o) 

The  divergence  of  this  part  is  evaluated  with  the  held  of  the  relation 
V-(a  b)  =  (V»a)b  +  (a-V)b: 


v.b}o)  „  ln(^)  .  3  jriH  [v  (*i!li)  ♦  lEfZl  v^o)]  -  3  b<°>b<°>.v  (JSilll  .      (49) 

B(o)2  2  g^(o)   -  -       -  B(o)2 


The  6-component    is 


Ve.v.ff(D    =   _3/r2   B(o).v    rj (o)j  (5Q) 

-     -  =]  -         -     B(o)2 


and  the  projection  along  the  magnetic  field  is 


B(o).^.ff(D  .  _2B(o).Vn(1)  +  3  _n^i  B(o).v[B(o)2j  _   (51 

B( 


"  -1  2  R(0)2  - 


The  part  of  l^1 '    proportional  to  uj1 ^  is  defined  by  ffp1 ^ 


41)  =  Ul(1)[W(0).B(°)  -  B<°>.W<°>) 


uj1)n(o)'[2V*(o)V<|,(o)    +    -2X(0)(V9B(0)    +   B(0)V6)    -   2rVo)2VeVe]    .       (52) 
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This  gives  for  the  divergence  of  the  stress  tensor 

V-B*1'  -  2  V^^V.(Ul(1)fi(0),V*(0))  ♦  u}1  V0)'V(|V^°>|2) 

♦  2Vry1(l)n(0),|V^(0)|2/r'  ♦  Vep} 1  ^(o)  '  x(o)rB(o) 

♦  lQrB{0)-l(v\])n{0)'rx{0))    ,  (53) 
with  the  projections  along  Ve 


ve.v.ff^)  =  -Lb^-vU^V^Vx^)        (5*0 

r 


and  along  B(o) 


B^-V-E^  m  V_L_1__   B(0).V(r2/|V0)|2)  +  X_!_  n^V^.yJr2^1*)  .   (55) 


The  remaining  part  of  E^1'  proportional  to  y},   is 

B(1)  .  uJ1)qCo)'[2x(o)2b(o)b(o)  .  r2B(o)2x(o)^(o)ye  +  y_eB(o))]  .    (56) 

The  divergence  yields 

7.*(1)  -  2ui1)n(0)'x(°)2[V(B(o)2/2)  *  J<°>xB<°>]  ♦  2b(°V^.V(u51)«(0),X(°)2) 

-  V6  B(°).V(ui1)"(0)VV0)2X(0))  *  2Vry51)«(0)'B(0)2X(0)2/r   (57) 


with  the  projections  along  Ve  and  3 
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(o) 


ye.v.iV)  -  ±-B{0)-ll2vi])n{0)\{0)3   -  Mii)aCo),P2B(o)2x(o)]  (58) 


and 


B(o).Zff(1)  „  b(o)2q(o)'x(o)2b(o).Zu(1)  _         (59) 


The  9-component  of  the  momentum  equation  then  gives  the  result 


-L  B(0).V(x(1)'r2n(0))  =  Ve-[j(o)*B(1)  +  J^xB^  ♦  V-B<1)]  , 
r2 


where  we  evaluate 


Ve.J(o)*B<1>  =  -<(°''/r¥0).y(1) 


and 


v9.j(1Wo)  -  JL  B(°).vx(1) 


Combining  this  with  the  results  of  eqs.   (50),  (5M)  and  (58),  we  obtain 
an  equation  of  the  form 


B(o)-V  {    }•=  0 


This  allows  us  to  solve  for  x 
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.0)  =  5(1)^(0)^  +  x(o)yn  +  x(n«fi(o)r2  +  3x(0)n(1)  „  ujl>B(o)'p2x(o> 


B(o)2 


-    2yPMo)Y°>3    +    y(1)„(0)'r2B(0)2x(0)  (6Q) 


with  a  free  surface   function   x^(<l^°^)« 

The   projection  of   the  momentum   balance  along  B^0'    is   of  the   form 

B(o).vr[-pO  Wo)2    -   2fl(0V1)p(0)/H    +    B^.VeA(1)yo).V(r2fl(o))    ♦    B(o).Vp(1) 


B(o).j(o)xB(1)    +   B(o).Vff(1) 


We  move  JxB^1^  to  the  left-hand  side  and  define  E  by  the  relation 


p(o)RT(o)  E  .  B(o).[(pv.v  V)(D  +  y  (1)-,  .  9P^1  B(o).v^(1)  p 


which  eventually  assumes  the  form 
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B(o).v  |p!H  +  rt^-RT^V"  .  diill  ,(D'  -^(0)'  ,.(D 


p 


RT(o)        RT(o)   c      p(o) 


~[^j,*(0)  H    l~p^     T<°>       4    H  rt<°> V 


R(0)   n   ™(1 ) 

T(1)_T(o)'  (1)    f.  (nS  ~  If 

I I * B(o).Vp(o)  =  — — — -  .  (61) 

T(o)p(o)    "    "        p(o)rT(o) 


Here  B(o).V.ff(1)  is  given  by  eqs.   (51), (55)  and  (59).   It  follows  that 
the  quantity 


r(°).\7.b>(1  ) 

^(1).t(o)'„,(1) 


T"  -r    *    B^-V  in   pCo)      (62) 


p(o)RT(o)        T(o) 


must  have  a  zero  flux  surface  average.   This  integrability  constraint 
is  used  to  determine  x'1   ,  which  occurs  linearly  in  Er°'«7'K|,   eq. 
(51),   through  n^1\  eq.   (48).  Note  that  the  dependence  between  A1 
and  i|/1^  is  linear.   The  explicit  equation  (61)  determines  then  p  ^   in 
the  form  of  a  modified  Bernoulli  law. 

The  first-'order  system  is  completed  by  the  Vy         component  of  the 
momentum  equation.  This  yields 
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V*(0).Vr[-p(1)rfl(0)2  -  2fi(o)v(1)p(o)/r]  +  V^o)-Vp( 


=  V^(o).[J(1)xB(o)  *  j(°^B(1)  +  V.B(1)]  .       (63) 


We  evaluate  the  Lorentz  force  terms 


V*(o).J(1)xB(o)  =  -x(0)/r2  V*(o).Vx(1)  -  |V*(o)|2/r2  A%(1 


V^°>.j(°W1>  m   .x(D/r2  |^(o)|2x(o)'  -  v*(1).V^0)/r2  A%(0)  . 


and  the  pressure  gradient 


V*(o).Vp(1)  =   V*(o).V(p(o)RT(1)  +  p(1)RT(o)] 


M  )    Ml    Ml'         ( 1  ) 
Inserting  the  results  for  "P   ,  x   »  *      and  Pi   we  obtain  an 

elliptic  Grad-SchlQter-Shafranov  type  equation  for  the  poloidal  flux 
(1)'. 


* 


V*^|2/r2  aV1)  =  G(r,A(1)',x(1)^(1),P(1)^(o),"0  ,   (64) 


where  G  contains  ijr  '  and  first  derivatives  of  ty  and  zerc-order 

quantities. 

Finally,   we  collect  the  results  for  the  first~order  corrections. 

1  "  M  )       "  M  ) 
There  are  three  free  surface  functions,  namely  <J>C,TV  '  and  x   •   We 

absorb  these  functions  in  the  zero«-order  equilibrium  quantities  and  can 

therefore  neglect  them  in  the  first  order.   Ohm's  law  provides  a 

relation  between  the  potential,  the  poloidal  flux  and  the  velocity  with 
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the  Hall  term  present.  The  energy  equation  yields  a  linear  relation 
between  "P  and  \\>  ,  and  \  *■  and  the  e-component  of  the  momentum 
equation  gives  a  similar  relation  for  the  toroidal  magnetic  field  x  • 
The  B^0'  component  of  momentum  balance  yields  a  solvability  constraint 
which  is  used  to  determine  for  Av  .  In  its  explicit  form  this 
equation  determines  the  density  p^  as  a  function  of  i|/  .  Finally, 
the  normal  component  determines  the  poloidal  flux  through  an  elliptic 
equation,  with  the  boundary  condition  that  i|/   vanish  at  the  surface. 
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V.   SECOND  ORDER 

In  second  order  in  e  we  admit  particle,  momentum,  energy  and 
current  sources.  We  also  allow  the  lowest-order  equilibrium  to  evolve 
on  the  diffusion  time  scale.  In  this  order  the  mass  and  energy 
diffusion  are  determined.  The  second-order  quantitites  are  not  that 
interesting  by  themselves;  they  are  therefore  only  evaluated  to  the 
extent  necessary  to  determine  diffusion.  They  averaging  procedure 
perfectly  describes  the  diffusion  and  makes  the  equations  simpler.  At 
this  point  the  time  evolution  of  the  free  surface  functions  has  to  be 
determined,  too,  in  order  to  close  the  system. 

To  maintain  steady  state,  there  must  be  a  mass  source  ps  in  the 
continuity  equation: 


V.(pv)<2)  =  p?  -±   p(Q)  .  (65' 

9t 


Rather  than  introduce  a  second-order  stream  function  for   (pv)1  '  we 
decompose 


(pv)<2)  =  -  "-*  °   ♦  VB<°>  -  p(0)v(2)  ♦  p"M"  ♦  p(2Mo)  ,   (66) 

lz*(0)l      -     -     -     - 

where  u  represents  the  normal,  diffusive  part  and  v  the  component  along 
the  field;  the  toroidal  component  is  not  essential  and  may  be  omitted. 
We  insert  the  representation  (66)  into  the  continuity  equation  and 
obtain 
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UV<J,(o) 


at       |y_^(0)|2 


This  allows  us  to  solve  for  v  up  to  an  arbitrary  function  of  i|r0'. 

The  condition  that  the  function  v  be  single-valued  yields  the 
integrability  constraint 


r  ji_  [S  _  ip^]  +  _d_  r  _d^     m 

J  B(o)  lP    at    d  (o)  J  B(o) 
p  p 


The  zero-order  profile  functions  p,  T^0',  ?r0'  and  x  depend  on  the 
flux  i|/o)  and  explicitly  on  the  time,  i.e.,  f  =  f  U(o)  (r,  z,t)  ;t ) .  The 
time  derivative  is  then  given  by  the  evolution  of  <|/0'  but  also  by  an 
explicit  dependence 


9f  .   3f   3ij/o)  +  3f , 
3t   3|j)(o)   3t     3tlex 


dil(o)    3.(o) 
The  first  part  is  no  longer  a  suface  quantity  since  — - depends 

0  t       0  t 


on  r  and  z,   but  the  explicit  part  — |ex  taken  at  ij/   kept  constant 


at  I 

still  is.   Instead  of  solving  very  complicated   integro-differential 


equations  for  the  temporal  evolution  of  the  four  profile  functions,  the 
explicit  time  derivative  are  expressed  in  terms  of  3ijr0  /3t  and 
equilibrium  quantities  of  known  order  and  the  evolution  of  ijr0  is 
determined  by  the  time  derivative  of  the  Grad-Schluter-Shafranov 
equation.   The  density  is  given  by  eq.   (33)  and  has  a  time  derivative 


-3U- 


3p 


(o) 


3t 


P(o)  [[p;  +  if  («!!!!] 

p    2   RT(o)  ,^o)         3t     9t  ^  /P    2  3t  I  RT  jext 


]  3*(o) 


-2  3  r^o)2 


The  surface  average  of  the  continuity  equation  then  yields 


,(o)  an(o) 


<p(o)>/;0  |PLv  +  <r2p(o)>  Q* 


3tlex 


RT 


(o)  at  lex 


2RT(o)2  at  lex 


<rVu'> 


i  ,  ■»   P  ,i,(o)   „2      0(o)2      ,  ai,,(o) 

4  (Jw)  fM>  *  <ps-p(0)  [-4—*  V  (2-ttt)  .,„,]  ^>  •   (68) 


,1» 


(oV 


RT(o)  ,^°>J  at 


In  (68)  J  is  the  Jacobian  of  the  transformation  from  r,9,z  to  ^,e,ep, 
flux  coordinates.  Next,  we  discuss  Ohm's  law,  which,  including  a 
source  term  Js,  reads 


.(1) 


n(2)j(o)  +  (n(2)   (2))j(o)  +  ("±_   JxB)(2)  +  j* 
1  "        II    1   "II       P   "- 


-1*{2)    +  "oEole  +I(0)x9(2)  +  v<1W1>  ♦  v<2>xB<°>  --^A(o)  ,   (69) 


where  r  E   is  the  constant  toroidal  field  at  r  =  rQ.  The  Hall  term 
assumes  the  form 


ni'V'l 


(n„/p  J.B)<2>  -  „(')/p(o)(j(o),bO)  .  j(D,b«")  -   "       j"»,B<°> 
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The  electric  and  magnetic  fields  are  expressed  by  scalar  and  vector 
potentials  where  A  has  the  form 

A  =  4)Ve  +  VAxVe 


with  the  gauge  V«A  =  0  .  In  order  to  obtain  a  representation  for  B 
consistent  with  that  of  eq.  (19),  for  a  given  x  the  function  A  has  to 
satisfy  A  A  =  ~x  with  the  boundary  condition  that  A  =  0  on  the  plasma 
surface.  We  consider  the  projection  of  Ohm's  law  along  B^0'  and,  using 
the  relations 


9A(0) 


B(o). «  J-  (V^0).VA(0J   ♦  X(0)  *(0J) 


r2 


b(o).,(1)xb(D  =  -b(D. (^1)^(0)  3P^  Zl|)(o)   +  ^(D) 

3ij/o) 


and 


b(D.v^D  =  -b<°).v  (a(o)V1)2  *  ♦(J)'  ♦">)  , 


we  obtain 
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.(o) 


n(2)[x(o)'B(o)2  +  x(o)  3Pi  +  ja.e(o)  .  .^(o)^  .  x(o)#(o)/rJ 

3r0' 


(71) 


-V^0).VAFt/r2  =  B(0).vU(2)-n(0U(2)-fi(0)^(1)2/2  ■ 


(1  )',,.(D 
c 


The  flux  surface  average  yields, using  the  relation 


,  LJL-*J°) 


■  Vi|ivu|.VA  t/r 

(o)  ~     -  >L 


d/0)  i 


dl   3X 


:o) 


B 


(o)   3t   ' 


and  after  differentiating  with  respect  to  ^ 


(o) 


3x 


3t   ex  3t 


(72) 


(o) 


♦  <1  {j[n(2>(x(o)  V°>2  ♦  X<0)  *£1L)    ♦  J8.B<°>  - 


,2E   (o)     /  n 


34; 


(o) 


♦ir"it(o)>  ■ 


which  states  that  we  can  express 


3X 


(o) 


3t 


3* 


(o) 


through  — and  equilibrium 

ex         ot 


quantities. 

Next  we  consider  the  9-component  of  Ohm's  law,  which  yields  after 
multiplication  by  r  q^0'   and  the  definition 


-37- 


Hl    .   p(o)jn(2)    3P(0)    1^(0)1: 


3lt)(o)       B(o)2 


*   n{f   j£L    [x(o)  V°)2    ♦    x(o)   IP^]    t   r2jS.7e   .  ,       (?3) 

II       B(o)2  3^,(0)  -     "  °   ° 


H}    ♦  B<0>.7{n<1)Vn   *   nin[-x(o)Vn   *   x(1)H   ^   P(o)  ^  -  u   .      (7H 

9 1 


This  equation  establishes  a  linear  relation  between  u  and  3ijr°v3t.  If 
3i|r°'/3t  is  given  then  y  is  known.  The  mass  diffusion  is  given  by  the 
surface  average 


<u> .  <Hl  ♦  p<°y°>> .  -  Uo)  m  (p^-^/M 

,c      J      J  b£o)        3t     J  B<o) 


relating  the  total  diffusion  to  the  sources  and  the  equilibrium.  With 
the  continuity  equation  the  mass  diffusion  through  a  surface  is 
determined  by  the  volume  integral  of  the  source  terms  ps  and  -3p^°'/3t. 
We  turn  to  the  energy  equation,  eq.  (3'')»  and  take  only  the 
surface  average.   The  left-hand  side  is  of  the  form 


LHS  =  [c  pv-VT  +  pV-v](2)  +  c  p(o)  II 


and,  furthermore, 


■38- 


LHS 


=    Cv[-T(0)'u    +    X(1)'b(0).Vt(D]    +    RT(0)[p?    - 


9tt 


(o) 


at 


RT(o)  v^(0).vP(0) 

7^~  ("M    \^{o)\2 


+    v    B(o).Vp(o)] 


-RT(0)^(1)»B(o).v(p(1)/p(o)}    .    RT(1>x(nV°>.Vp<°>    ♦    cvp(0)    ill 


(o) 


With  the  relations 


,(o)        " 


(o)2„2 


-   rt(0)    I  J!_  v  B<°>.V   £n  p<°>   =   |  -?*      °  •      B(Q).7v 


,(0)  2 


and 


d!.  2n 
r       V  • 


uV* 


:o 


B(o) 


,(o)|2        aih(o)    J 


dX,       2 


v   -    2 


|V^0J|^        9<|/o;    J    B^ 


di 


u  rij; 


(o) 


(o)   ■  'I  B(o)    |^(o)|2 


and   the  definition 


H,    .    <CvT(0)'(p(0)    2^1  .    u)    .    „<•>    J    (Jm)^(o)    +    „(0)    ^  -!^J! 


♦  R  ^D'o^b^-VT^ 


(75) 
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we  obtain  for  the  surface  average 


<LHS>  =  c. 


3T 


<p(°»>  +  5 <r2B(o)-Vv>  +  H?  .    (76) 

ex  2     - 


3t 


We  now  consider  the  right-hand  side  of  the  energy  equation: 


RHS 


.  -y.^(2)  +  (ff:v  v)(2)  +  [J'(E+VxB)](2)  +  E   , 


where 


ff:7  v 


(2>  =  1t~  m){2) 

2 


!°!l  Tr  W<°>2  *   I!L-  Tr  (B<°V°>)2  -  [^l!2.  (B.W.B)2](2' 


2B< 


The  first  term  has  with  eq.   (28)  the  form 


R,  =  p(2)r2(^0),)2|V^°)|2 


and  the  second 


4zha{0)')z\l*l0)\* 


The  third  term 


u(0) 

R,  =  -i  — I (B.W«B)(1)2 

3     2  2B(o)2 


is  determined  from  the  first-order  calculation,  eq.    (48).   Then  the 
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viscous  energy  dissipation  is  known.   The  heat  flax  contains  in  second 
order  three  terms 


V2)  -  (K  Vu  T)(2)  ♦  K<2M°>  V0)  -  (K/x  15!)  (2)  . 

ir  ii         i       _         A  B 


For  the  heat  flux  we  separate  the  contributions  along  B^0'  and  ET   a 
obtain 


,K/\   X 
■  -   B(o) 


-  V-2(2)  =  b(°)-VQ0  +  B^-VQ,  +  V.(k|2Mo)  V0))  +  1.{J1-1 Vt(1}xV( 


Since  the  term  B^u;«VQ   does  not  contribute  to  the  surface  integral 


,<o>., 

only  Q1  is  required 


(1)„(o) 


Q,    =   -    K/X      X  T<°>'    ♦    T<°V1>V°>    .  (77) 


B(o) 


If  the  heat  generation  term  Q   is  included   it  contributes  to  the 
surface  integral  as 


e  ^(o)  B(o)2  -    "L  2  ^2me;  B(o) 


We  define 


•41 


H3  =  ^.VQ,  ♦  1  (Jk|2¥o)'|V^o)|2 


1 


»¥ 


J  o(0)  -        '   J.(1  ) 


B( 


.* 


+  R,  +  R2  +  R3  +  n(2)J(0)2 


n]fj];)2  +  J-J3  ♦  Eg  ♦  Qe> 


and  obtain  from  energy  balance 


(oK  3T(o)    +a^f<r2B(o).Vv>  =  H 

ex     2      ~    -      3    2 


cv  <p^'> 
v         3t 


(78 


Finally,  the  momentum  balance  is  considered.  It  is  sufficient  to 
evaluate  the  surface  average  of  the  9-component.  Then  no  second^order 
fields  are  involved.   The  LHS  contributes 


LHS  =  p(o)  12_L  +  V6-(pv.V  v)(2) 
3 1 


-J-  [p(o>  1^+  (pI)(2),v(»'  ♦  (pv)C).Vv(')] 
p2         3t 


After  multiplication  by  r  and  taking  the  surface  average  we  define 


,(o) 


,(o) 


at  |v^(o)|2 


>     (79) 


The  RHS  has  as  one  term 
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ve.(W2)  =  JL  [b(°).7(x(2)-x(0)^(2))  ♦  B^.7X(n] 


The  pressure  tensor  contains  three  terms 


.(2)    .   u(2)w(o)    +    u(2)B(o).w(o).B(o) 


K^2)=    uj1)(W.B-B-W)(1  }    +    uin(B   B.W.B-B-W.B   B)(1) 


and 


g(o)g(o)  g   g 

,(2)   -   40)[Ul   "   3  ~       ~        ]    (B-W.B)(2)    -    3    (  — )(1)    (B-W.B)(1)] 

;3  *  B(o)2  ~  -  ~  B2  ~  -  ~ 


m   ,  R(0)R(0))(B-w.B)(1) 


The  only  contribution  from  ffl2^    is,  with  eq.      (48), 

<r2ve-v.^2)>  =  -3<B(1).V(x(o)n(1)/B(o)2)> 

containing  only  known  quantities.  The  remaining  part  of  the  pressure 
tensor  involves  only  zero*-  and  first  order  quantitites;  its  evaluation 
is  straightforward.   With  the  definition 


<-U3- 


H5  =  <3(1)-VX(1)  +  r2Ve-V(ff{2)+  1{22)) 


-   3  B(1  ^VCx^V1  ]/B{o)2)   +  r2Ve-Es>  (80) 


the  averaged  momentum  balance  yields  the  result 


9"  °  I   <r2p(o)>  -  n(0)<r2B(0).Vv>  =  H,   -   Hn        (81) 
3t   ex 


Let  us  first  discuss  the  stead"y  state  solution.  Ohm's  law,  eq. 
(72),  relates  the  source  to  the  fluxes;  eq.  (74)  determines  the 
perpendicular  mass  flow  in  dependence  of  the  source  and  the 
equilibrium:  The  energy  outflow  through  a  surface  is  also  determined  by 
the  sources  and  the  equilibrium,  eqs.   (78)  and  (81). 

The  time  evolution  of  the  profiles  is  given  by  their  explicit  time 
dependence  and  the  evolution  of  the  flux  and  is  computed  iteratively. 
The  explicit  time  dependence  is  assumed  to  be  given.  Then  the  final 
equation  needed  to  close  the  system  is  the  time  derivative  of  the 
Grad-SchlUter«iShafranov  equation : 


A*  1*^-  n  ±    (.2  9p!fl  +  „(0)(0)' 

at    3t  ^"  3l|)(o)   * 


-  -     3     fr>2    9P(0)    +    y(o)    (o)")        _   9*^1  _9_   rp2    9p^  +      (o)    (0)'l      (g2) 
3t   l"      9^(o)        X        X  Jex  3t      ^(o)    l        3„,<o) 
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The  RHS  is  linear  in  3<Jj   /3t.  The  boundary  condition  requires  that 
the  tangential  electric  field  to  be  constant  at  the  plasma  surface 


3*(0) 

— =  r^En  =  ct  .   at  the  suface. 

3t     °  ° 


Since  eq.  (82)  is  a  linear  generalized  differential  equation, 
according  to  Ref.  [131.  it  requires  an  additional  boundary  condition  at 
the  origin  3<Jr°V3t  =  c  at  the  origin. 

Next  assume  3iK°V3t  to  be  given.   Ohm's  law  establishes  a  linear 
relation  between  u  and  34>(o)/3t,  eq.   (7*0,  x(o).  eq.    (72).   By  the 

s   a0(o) 

continuity  equation  <p  -  — >  and  momentum  balance,  eqs.   (78)  and 

gT(o)                       afl(°^ 
(81),  1    is  related  linearly  to  LY.   From  the  continuity 

at  i ex  at  i e* 

equation  eq.   (68),  and  either  energy  or  momentum  balance  the  remaining 

3p1,      .  _    3T(o 
_|ex  arld,  say,  _ 


3o1             3T^0^ 
two  profiles  \--—  Lv  and,  say,  )LV  are  determined.   Eventually 

-  | ex  gt    | ex 

the  explicit  profile  change  is  expressed  through  34r°V3t  and  known 
equilibrium  quantities. 

The  profiles  f  =  f(ijj'   ,t)  are  represented  by  cubic  splines,  where 

the  expansion  parameters  depend  on  time 

f  =  I   an(t)en(*(o))  . 
n 


This  leads  to  systems  of  ordinary  differential  equations  for  the 
coefficients,  which  can  be  solved  using  standard  techniques. 
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VI.   DISCUSSION 

The  general  compressible,  resistive  and  viscous  MHD  equations 
serve  as  a  basis  for  the  simulation  of  tokamak  evolution  on  the 
diffusion  time  scale.  Classical-type  transport  coefficients  are 
employed  in  this  single-fluid  model.  A  strong  magnetic  field  requires 
a  more  general  stress  tensor  than  that  from  ordinary  fluid  theory. 
This  general  tensor  can  be  derived  with  only  covariance  and  symmetry 
properties  being  used  (Ref.  [2]),  implying  that  the  Braginskii  form 
derived  from  microscopic  theory  with  genei"al  viscous  coefficients  is 
the  most  generally  admissible  stress  tensor.  The  transport 
coefficients  are  grossly  disparate  in  magnitude.  This  fact  prevents 
the  direct  numerical  solution  for  the  very  small  perpendicular 
diffusion  superimposed  on  large  flow  within  a  magnetic  surface. 

The  transport  coefficients  parallel  and  perependicular  to  the 
magnetic  field  differ  by  the  factor  u.t^  and  a  systematic  expansion  is 
therefore  applied,  using  e  =  (oj^t^)  as,  expansion  parameter.  The 
steady-state  solution  yields  in  lowest  order  a  Grad~Schluter*-Shafranov 
equation  allowing  a  purely  toroidal  flow  depending  on  the  flux  alone, 
and  a  temperature  constant  on  a  flux  surface.  In  first  order  i. 
poloidal  flow  becomes  necessary,  but  this  flow  is  within  a  magnetic 
surface.  In  second  order  in  e,  this  being  the  time  scale  for  the  time 
evolution  of  the  profiles  perpendicular  mass  and  energy  transport 
occur.  Introducing  surface  averaging  then  requires  the  solution  of 
general  Grad-Schluter*<Shafranov  equations  together  with  the  time 
integration  for  the  zero^order  flux,  yielding  a  11/2  D  algorithm 
alternating  between  2D  equlibria  and  1D  transport.  A  numerically 
tractable  scheme  is  thus  obtained  which  can  be  solved  efficiently  by 
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means  of  the  finite-element  method.  This  self-consistent  model  is 
Suitable  for  realistic  description  of  experimentally  relevant 
configurations;  (1)  large  toroidal  flow,  (2)  small  poloidal  flow  and 
poloidal  variation  on  a  flux  surface  for  the  quantities,  and 
(3)  enhanced  perpendicular  energy  and  mass  outflow.   Especially  the 

gyroviscous  part  of  the  pressure  tensor  associated  with  the  toroidal 

( n) ' 
velocity  variation  ftv     produces  such  an  enhanced  loss.   General 

boundary  conditions  and  source  terms  can  be  applied  as  discussed  in 

detail  in  the  paper. 

Our  model  applies  strictly  for  collision-dominated  plasmas,  and  so 

the  1  1 /2  D  scheme  introduced  is  appropriate.   For  very  hot  plasmas, 

especially  for  tokamaks  with  anomalous  transport,  the  diffusion  rate 

and  profile  evolution  time  scale  of  the  order  of  (id^t^)   z   10'     are 

too  small  compared  with  experiments.   It  may  be  argued  that  the  set  of 

macroscopic  equations  is  still  appropriate;  however,   the  transport 

coefficients   need   to   be   adjusted  with  a  "cutoff"  value  for 

(wjTj)   z   10"3  and   enhanced   transport   coefficients.    Such   a 

phenomenological  approach  emerging  from  a  rigorous  treatment  of  the  set 

of  dissipative  equations  should  complement  the  more  heuristic  1D  and  1 

1/2  D  transport  simulations.   It  seems  worthwhile  fully  to  explore  the 

solutions  of  such  a  "generalized  classical  model"  instead  of  purely 

adjusting  transport  coefficients  in  a  non'-self'-consistent  model. 
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